#!/usr/bin/env python3
# -*- coding: utf-8 -*-

__author__ = 'Nan Wang'

################常数或建议数值
f_s = 0.05 #溶解性的惰性COD占总COD比例,0.05~0.1,市政污水建议去0.05
f_A = 0.3 #颗粒性惰性组分比例,0.2~0.35,对于城市污水取0.3
f_COD = 0.2 #易降解COD比例,0.15~0.25
f_B = 0.3 #进水可过滤无机物质,进厂污水取0.3,初沉池污水取0.2
S_orgN_AN = 2.0 #mg/L,出水有机氮
Y_COD_abb = 0.67 #可降解COD产泥系数,降解每gCOD形成的生物量COD
b = 0.17 #d^-1,15℃衰减系数
SVI = 120.0 #120L/kg,100-150,污泥体积指数
t_E = 2.0 #h,设计浓缩时间不能太长，防止污泥溶解和二沉池反硝化使沉淀污泥悬浮
R_m2a = 4.00 #膜池2好氧池回流比
R_D3an = 2.00 #缺氧池2厌氧池回流比
C_0 = 2.00 #mg/L,混合液剩余DO值
C_S = 9.17 #mg/L,标准条件下清水中饱和溶解氧
T_TS = 20 #摄氏度，计算混合液温度(标况)
T_summer = 35 #摄氏度，夏季温度
h_p = 0.4 #m,管道阻力
h_A = 0.4 #m,曝气器水头损失
h_delta = 0.02 #m,每升高1℃需补偿压力值
S_NO3_ZB = 0.0 #mg/L，设定进水硝酸盐氮为0
miu_A_max = 0.47 #d^-1,最大比生长速率


###一 设计流量######
Q_d_Knoz = float(input("日平均进水水量(m3/d)"))
Q_h_Knoz = Q_d_Knoz / 24
if Q_h_Knoz <= 13:
    Kz = 2.7
elif Q_h_Knoz >= 2600:
    Kz = 1.5
else:
    Kz = 3.5778 * Q_h_Knoz**(-0.112) #变化系数
Q_d_max = Q_d_Knoz * Kz
Q_h_max = Q_d_max / 24
print("变化系数 %.2f\n日最高进水水量 %d m3/d" %(Kz,Q_d_max))

###二 进出水水质及平衡#######
#2.1输入参数
print('------输入进水水质------')
C_COD_ZB = float(input('进水化学需氧量(mg/L)'))
C_BOD5_ZB = float(input('进水生物需氧量(mg/L)'))
C_P_ZB = float(input('进水总磷(mg/L)'))
C_TN_ZB = float(input('进水总氮(mg/L)'))
C_SS_ZB = float(input('进水悬浮固体(mg/L)'))
T_C = float(input('设计温度(℃)'))
print('------输入出水水质------')
S_COD_AN = float(input('出水化学需氧量(mg/L)'))
S_BOD5_AN = float(input('出水生物需氧量(mg/L)'))
S_TP_AN = float(input('出水总磷(mg/L)'))
S_TN_AN = float(input('出水总氮(mg/L)'))
S_NH4_AN = float(input('出水氨氮(mg/L)'))
S_SS_AN = float(input('出水悬浮固体(mg/L)'))

#2.2碳平衡
print('------碳平衡------')
X_TS_ZB = C_SS_ZB #进水可过滤物质
X_COD_ZB = X_TS_ZB * 1.6 * (1-f_B) #颗粒性COD(可过滤物质COD),有机干物质颗粒按1.6gCOD/oTS计
S_COD_ZB = C_COD_ZB - X_COD_ZB #可溶解性COD
S_COD_inert_ZB = f_s * C_COD_ZB #溶解性惰性组分
X_COD_inert_ZB = f_A * X_COD_ZB #颗粒性惰性组分
C_COD_abb_ZB = C_COD_ZB - S_COD_inert_ZB - X_COD_inert_ZB #可降解COD
C_COD_la_ZB = f_COD * C_COD_abb_ZB #易降解COD

X_anorg_TS_ZB = f_B * X_TS_ZB #进水可过滤无机物质(仅算数,进水颗粒性COD没有直接用)

print('可降解COD %.2f(mg/L)\n易降解COD %.2f(mg/L)'%(C_COD_abb_ZB,C_COD_la_ZB))

print('进水化学需氧量 %.2fmg/L=可降解COD %.2fmg/L+溶解性惰性组分%.2fmg/L+颗粒性惰性组分%.2f\
mg/L'%(C_COD_ZB,C_COD_abb_ZB,S_COD_inert_ZB,X_COD_inert_ZB))

#2.3出水氮平衡
print('------氮平衡------')
S_TKN_AN = S_NH4_AN + S_orgN_AN  #mg/L,出水凯氏氮
S_anorgN_UW = S_TN_AN - S_TKN_AN #mg/L,出水硝酸盐氮

###三 硝化菌泥龄#######
B_d_COD_Z = Q_d_Knoz * C_COD_ZB / 1000 #mg/L,COD日负荷
def PF_F(B_f):
    if B_f <= 2400:
        PF_f = 2.1
    elif B_f > 12000:
        PF_f = 1.5
    else:
        PF_f = 2.1 - (B_f - 2400) * 0.6 /9600
    return PF_f
PF = PF_F(B_d_COD_Z) #硝化反应系数
t_TS_aerob_Bem = PF * 1.6 / miu_A_max * 1.103**(15-T_C) #d,硝化菌污泥龄
print("硝化反应系数 %.2f\n硝化菌泥龄%.2fd"%(PF,t_TS_aerob_Bem) )

###四 反硝化体积比例VD/VBB#######
#V_D_over_V_BB = 0.2 #前缺氧池/曝气池体积比例
#C_COD_dos = 0 #mg/L,外加碳源化学需氧量

COD_dos_name = input("请输入甲醇,乙醇或醋酸之一") #投加碳源类型
if COD_dos_name == '甲醇': 
    Y_COD_dos = 0.45
elif COD_dos_name == '乙醇' or '醋酸': 
    Y_COD_dos = 0.42
F_T = 1.072**(T_C-15) #内源呼吸的衰减系数


def V_D_over_V_BB_F():
    C_COD_dos_f = 0     #外加碳源化学需氧量
    V_D_over_V_BB_f = 0.2
    x_f = 0 
    while x_f < 1:
        #4.1污泥产量的计算
        t_TS_Bem_f = t_TS_aerob_Bem / (1-V_D_over_V_BB_f) #设计污泥泥龄
        X_COD_BM_f = (C_COD_abb_ZB * Y_COD_abb + C_COD_dos_f * \
            Y_COD_dos)/(1 + b * t_TS_Bem_f * F_T)     #生物体中的COD
        X_COD_inert_BM_f = 0.2 * X_COD_BM_f * t_TS_Bem_f * b * F_T  #剩余惰性固体
        US_d_C_f = Q_d_Knoz * (X_COD_inert_ZB / 1.33 +(X_COD_BM_f + \
            X_COD_inert_ZB) / (0.93 * 1.42) + f_B * X_TS_ZB) / 1000 #污泥产量
        #4.2反硝化硝态氮浓度计算
        S_NO3_AN_f = 0.7 * S_anorgN_UW #出水硝态氮
        X_orngN_BM_f = 0.07 * X_COD_BM_f #形成活性污泥的氮
        X_orgN_inert_f = 0.03 * (X_COD_inert_BM_f + X_COD_inert_ZB) #与惰性颗粒结合的氮
        S_NO3_D_f = C_TN_ZB - S_NO3_AN_f - S_orgN_AN - S_NH4_AN - \
            X_orngN_BM_f - X_orgN_inert_f #每日平均反硝化的硝态氮浓度
        #4.3碳降解的需氧量
        OV_C_f = C_COD_abb_ZB + C_COD_dos_f - X_COD_BM_f - \
            X_COD_inert_BM_f #碳降解的总需氧量
        OV_C_la_vorg_f = f_COD * C_COD_abb_ZB * (1-Y_COD_abb) +\
            C_COD_dos_f*(1-Y_COD_dos) #反硝化区易降解及外加碳源需氧量
        OV_C_D_f = 0.75 * (OV_C_la_vorg_f + (OV_C_f - OV_C_la_vorg_f) *\
            V_D_over_V_BB_f**0.68) #反硝化区总需氧量
        #4.4耗氧量和供氧量平衡
        x_f = OV_C_D_f / 2.86 / S_NO3_D_f
        if V_D_over_V_BB_f < 0.6 and x_f < 1:
            V_D_over_V_BB_f += 0.01
        elif V_D_over_V_BB_f >= 0.6 and x_f < 1:
            V_D_over_V_BB_f = 0.6
            C_COD_dos_f += 0.01

    return(C_COD_dos_f, V_D_over_V_BB_f, x_f, t_TS_Bem_f, X_COD_BM_f,\
           X_COD_inert_BM_f, US_d_C_f, S_NO3_AN_f,X_orngN_BM_f, X_orgN_inert_f,\
            S_NO3_D_f, OV_C_f,OV_C_la_vorg_f,OV_C_D_f)

result = V_D_over_V_BB_F()
C_COD_dos = result[0]   #mg/L，外加碳源化学需氧量
V_D_over_V_BB = result[1]   #缺氧池/曝气池体积比例
x = result[2]   #耗氧量和供氧量平衡
t_TS_Bem = result[3]    #设计污泥泥龄
X_COD_BM = result[4]    #生物体中的COD
X_COD_inert_BM = result[5]  #剩余惰性固体
US_d_C = result[6]  #污泥产量
S_NO3_AN = result[7]    #出水硝态氮
X_orngN_BM = result[8]  #形成活性污泥的氮
X_orgN_inert = result[9]    #与惰性颗粒结合的氮
S_NO3_D = result[10]    #每日平均反硝化的硝态氮浓度
OV_C = result[11]   #碳降解的总需氧量
OV_C_la_vorg = result[12]   #反硝化区易降解及外加碳源需氧量
OV_C_D = result[13]     #反硝化总需氧量

T_TS_D_Bem = t_TS_Bem - t_TS_aerob_Bem #d,反硝化菌泥龄

print('缺氧池/曝气池体积比例%.2f'%V_D_over_V_BB)
print('外加碳源化学需氧量%.2fmg/L'%C_COD_dos)
print('耗氧量和供氧量平衡%.0f'%x)
print('设计污泥泥龄%.2fd'%t_TS_Bem) 
print('污泥产量%.2fkg/d'%US_d_C)
print('每日平均反硝化的硝态氮浓度%.2fmg/L'%S_NO3_D)
print('碳降解的总需氧量%.0fmg/L'%OV_C)
print('碳反硝化总需氧量%.2fmg/L'%OV_C_D)
###五 除磷######
#5.1 生物处理与化学除磷量
C_P_AN = 0.7 * S_TP_AN  #mg/L,出水浓度
X_P_BM = 0.005 * C_COD_ZB   #形mg/L,成活性污泥的氮
X_P_BioP = 0.006 * C_COD_ZB    #mg/L,生物法除磷量
X_P_Fall = C_P_ZB - C_P_AN - X_P_BM - X_P_BioP  #mg/L,需要沉析的磷酸盐
Me_3plus = 1.5 * X_P_Fall / 31   #mol/L,化学除磷药剂投加量
#5.2 除磷污泥产量
P_dos_name = input('化学除磷药剂的选择,输入铁盐或铝盐')
if P_dos_name == '铝盐':
    X_P_Fall_Fe = 0  #折合铁盐投加量
    X_P_Fall_Al = 27 * Me_3plus     #折合铝盐投加量
elif P_dos_name == '铁盐':
    X_P_Fall_Fe = 55.8 * Me_3plus   #mg/L,折合铁盐投加量
    X_P_Fall_Al = 0   #mg/L,折合铝盐投加量   
US_d_P = Q_d_Knoz * (3 * X_P_BioP + 6.8 * X_P_Fall_Fe + 5.3 * X_P_Fall_Al) / 1000   #化学除磷产泥量
#5.3 污泥产量
US_d_r = US_d_C + US_d_P    #kg/d,剩余污泥量
M_TS_BB = t_TS_Bem * US_d_r   #kg,生物段保持的污泥质量
M_TS_D = V_D_over_V_BB * M_TS_BB    #kg,缺氧池污泥量
M_TS_aero = M_TS_BB - M_TS_D    #kg,好氧池污泥量
K_de = Q_d_Knoz * S_NO3_D / M_TS_D / 1000   #kgN/kgSS·d,反硝化速率
L_C = (C_COD_ZB - S_COD_AN) * Q_d_Knoz / 1000 / M_TS_aero #kgCOD/kgSS·d,好氧池COD负荷
L_B = (C_BOD5_ZB - S_BOD5_AN) * Q_d_Knoz / 1000 / M_TS_aero #kgBOD/kgSS·d,好氧池BOD负荷

print('剩余污泥量 %.2fkg/d\n反硝化速率 %.2fkgN/kgSS·d\nBOD负荷 %.2fkgBOD/kgSS·d'%(US_d_r, K_de, L_B))

###六 污泥浓度######
TS_BB = float(input('生物池污泥浓度(g/L)'))  #生物池污泥浓度

###七 生物池容积######
#7.1生物池容积
V_BB = M_TS_BB / TS_BB  #m3,曝气池的容积
V_an = 1 * Q_h_Knoz     #m3,厌氧池容积
V_D = V_BB * V_D_over_V_BB  #m3,缺氧池容积
V_aero = V_BB - V_an    #m3,好氧池容积
V_bioT = V_BB + V_an    #m3,总容积
HRT_an = V_an / Q_h_Knoz     #h,厌氧池水力停留时间
HRT_D = V_D / Q_h_Knoz  #h,缺氧池水力停留时间
HRT_aero = V_aero / Q_h_Knoz   #h,好氧池水力停留时间
HRT_bioT = V_bioT / Q_h_Knoz    #h,总水力停留时间
#7.2 回流比
RF = S_NO3_D / S_NO3_AN     #反硝化所需的回流比
RZ = RF - 1  #反硝化所需的内回流比
eta_0 = 1 - 1 / (1 + RF)        #反硝化最大效率
print('生物池总容积 %.2fm3|生物池总停留时间 %.2fh'%(V_bioT, HRT_bioT))
print('厌氧池总容积 %.2fm3|厌氧池总停留时间 %.2fh'%(V_an, HRT_an))
print('缺氧池总容积 %.2fm3|缺氧池总停留时间 %.2fh'%(V_D, HRT_D))
print('好氧池总容积 %.2fm3|好氧池总停留时间 %.2fh'%(V_aero, HRT_aero))
print('内回流比%.2f '%RF)
print('反硝化效率%.2f'%eta_0)

###八 需氧量######
#8.1 耗氧量物料平衡
OV_d_C = Q_d_Knoz * OV_C / 1000     #kgO2/d, 碳去除的耗氧量
OV_d_N = Q_d_Knoz * 4.3 * (S_NO3_D - S_NO3_ZB + S_NO3_AN) / 1000   #kgO2/d, 反硝化回收供氧量
OV_d_D = Q_d_Knoz * 2.86 * S_NO3_D / 1000   #kgO2/d, 反硝化回收供氧量
OV_h_aM = ((OV_d_C - OV_d_D) + OV_d_N) / 24     #kgO2/h，平均耗氧量
OV_h_max = Kz * OV_h_aM     #kgO2/h，最高耗氧量
print('碳去除的耗氧量%.2fkgO2/d\n反硝化回收供氧量%.2fkgO2/d\n反硝化回收供氧量%.2fgO2/d\n平均耗氧量%.2fkgO2/h'%(OV_d_C, OV_d_N, OV_d_D, OV_h_aM))
# 8.2 标准传氧速率
alfa = float(input('α(0.8-0.85)'))  #0.8-0.85,混合液KLa/清水KLa
beta = float(input('β(0.9-0.97)')) #0.9~0.97，混合液饱和溶解氧/清水饱和溶解氧
h_TB2A = float(input('曝气装置与池底距离(m)'))
h_tk = float(input('设计水深(m)'))
h_El = float(input('当地海拔高度(m)'))
E_A = float(input('氧利用率'))
O_t = 21 * (1 - E_A) / (79 + 21 * (1 - E_A))
P_a = (101325 - h_El /12 / 133) / 1000000 #Mpa,当地大气压力
P_b = P_a + (h_tk - h_TB2A) * 9.81  / 1000 #Mpa曝气装置处绝对压力
C_SW = 8.24 * P_a / 0.101325    #mg/L,清水表面饱和溶解氧
C_SM = C_SW * (O_t / 42 + P_b / (2 * P_a)) #mg/L,水下深度到池面清水平均溶氧值
FCF = alfa * (beta * C_SM - C_0) / C_S #AOR与SOR转换系数
SOR = OV_h_aM  / FCF #kgO2/h,标准传氧速率SOR
G_S = SOR / (0.28 * E_A)
V_GS_over_V_knoz = G_S / Q_h_Knoz
print('标准传氧速率SOR%.2fkgO2/h\n标准状况供空气体积%.2fm3/h\n气水比%.2f'%(SOR,G_S,V_GS_over_V_knoz))